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Abstract —We treat the problem of estimation of orientation 
parameters whose values are invariant to transformations from 
a spherical symmetry group. Previous work has shown that 
any such group-invariant distribution must satisfy a restricted 
Unite mixture representation, which allows the orientation pa¬ 
rameter to be estimated using an Expectation Maximization 
(EM) maximum likelihood (ML) estimation algorithm. In this 
paper, we introduce two parametric models for this spherical 
symmetry group estimation problem: 1) the hyperbolic Von Mises 
Fisher (VMF) mixture distribution and 2) the Watson mixture 
distribution. We also introduce a new EM-ML algorithm for 
clustering samples that come from mixtures of group-invariant 
distributions with different parameters. We apply the models 
to the problem of mean crystal orientation estimation under the 
spherically symmetric group associated with the crystal form, e.g., 
cubic or octahedral or hexahedral. Simulations and experiments 
establish the advantages of the extended EM-VMF and EM- 
Watson estimators for data acquired by Electron Backscatter 
Diffraction (EBSD) microscopy of a polycrystalline Nickel alloy 
sample. 

I. Introduction 

This paper considers estimation of parameters of distribu¬ 
tions whose domain is a particular non-Euclidean geometry: 
a topological space divided into M equivalence classes by 
actions of a finite spherical symmetry group. A well known 
example of a finite spherical symmetry group is the point 
group in 3 dimensions describing the soccer ball, or football, 
with truncated icosahedral symmetry that also corresponds 
to the symmetry of the Carbon-60 molecule. This paper 
formulates a general approach to parameter estimation in 
distributions defined over such domains. We use a restricted 
finite mixture representation introduced in jl] for probability 
distributions that are invariant to actions of any topological 
group. This representation has the property that the number 
of mixture components is equal to the order of the group, the 
distributions in the mixture are all parameterized by the same 
parameters, and the mixture coefficients are all equal. This 
is practically significant since many reliable algorithms have 
been developed for parameter estimation when samples come 
from finite mixture distributions 0, 0- 

We apply the representation to an important problem in 
materials science: analysis of mean orientation in polycrys¬ 
tals. Crystal orientation characterizes properties of materials 
including electrical conductivity and thermal conductivity. 


Polycrystalline materials are composed of grains of varying 
size and orientation, where each grain contains crystal forms 
with similar orientations. The quality of the material is mainly 
determined by the grain stmcture i.e. the arrangement of the 
grains, their orientations, as well as the distribution of the 
precipitates. Thus accurate estimation of crystal orientation of 
the grains is useful for predicting how materials fail and what 
modes of failure are more likely to occur Q. 

The mean orientation of the grain, characterized for example 
by its Euler angles, can only be specified modulo a set of 
angular rotations determined by the symmetry group associ¬ 
ated with the specific type of crystal, e.g. hexagonal, cubic. 
This multiplicity of equivalent Euler angles complicates the 
development of reliable mean orientation estimators. The prob¬ 
lem becomes even harder when the orientations are sampled 
from a region encompassing more than one grain such that 
the orientations cluster over different mean directions. In such 
a case, we would like to identify whether the orientations are 
multi-modally distributed and also estimate the mean direction 
for each cluster. 

In our previous work |[T}, we introduced the finite mixture of 
Von Mises-Fisher (VMF) distribution for observations that are 
invariant to actions of a spherical symmetry group. We applied 
the expectation maximization (EM) maximum likelihood (ML) 
algorithm, called EM-VMF, to estimate the group-invariant 
parameters of this distribution. In this paper, we develop 
a hyperbolic representation simplification of the EM-VMF 
algorithm that reduces the computation time by a factor of 
2. We also introduce a new group invariant distribution for 
spherical symmetry groups, called the (/-invariant Watson 
distribution, which like VMF is a density parameterized by 
location (angle mean) and scale (angle concentration) over 
the //-dimensional sphere. An EM algorithm is presented for 
estimation of the parameters, called the EM-Watson algorithm. 
Furthermore, mixture-of-C/-invariant Watson (mGIW) and von 
Mises-Fisher (mGIV) distributions are introduced to perform 
clustering on the (/-invariant sphere. An EM algorithm is 
presented for estimation of the parameters of the mGIW 
and mGIV distributions. We illustrate how the Generalized 
Likelihood Ratio Test (GLRT) can be used to detect the 
presence of multiple modes in a sample and how it can 
be combined with the EM algorithm for mGIW and mGIV 


distributions to cluster multiple orientations on the sphere. 

The performance of the proposed EM orientation estimators 
is evaluated by simulation and compared to other estimators. 
The EM orientation estimators are then illustrated on Electron 
Backscatter Diffraction EBSD data collected from a Nickel 
alloy whose crystal form induces the m3m [5] cubic point 
symmetry group. We establish that the EM orientation esti¬ 
mators result in significantly improved estimates of the mean 
direction in addition to providing an accurate estimate of 
concentration about the mean. Furthermore, with the extended 
mixture models, we are able to identify and cluster multi- 
modally distributed samples more accurately than the K-means 
algorithm. 

The paper is organized as follows. Section [D] describes 
group invariant random variables and gives the mixture repre¬ 
sentation for their densities. Section [III] specializes to random 
variables invariant relative to actions of the spherical symmetry 
group and develops the Q -invariant VMF and Watson distribu¬ 
tions along with EM-ML parameter estimator. The clustering 
methods based on the Cl-invariant distributions along with 
the GLRT are elaborated in Section IV The crystallography 
application and data simulation are presented in Section [V| and 
the experiment results are shown in Section [Vi] Section |VII| 
has concluding remarks. 


II. Group-invariant random variables 


Consider a finite topological group Q = {G\ 1 ..., Gm} of 
M distinct actions on a topological space X, G, : X -A X 
and a binary operation defining the action composition 
Gi * Gj, denoted G t G :] . Q has the properties that composition 
of multiple actions is associative, for every action there exists 
an inverse action, and there exists an identity action [6). A real 
valued function /(x) on X is said to be invariant under Q if: 
/(Gx) = /(x) for G £ Q. Let X be a random variable defined 
on X. We have the following theorem for the probability 
density /(x) of X. 

Theorem 2.1: The density function / : X — > M is invariant 
under Q if and only if 


3 h : X —> R. s.t. 

1 M 


/( x ) = 


(i) 


This theorem is a slight generalization of | Ij Thm. 2.1] in that 
the density h(.) is not necessarily the same as /(.). The proof 
is analogous to that of JT, Thm. 2.1]. 

Theorem 2.1 says that any density /(x) that is invariant 


under group Q can be represented as a finite mixture of a 
function and its translates /i(Gjx) under the group’s actions 
Gi £ Q. As pointed out in [Tj, Thm. 2.1 has important 
implications on (/-invariant density estimation and parameter 
estimation. In particular it can be used to construct maximum 
likelihood estimators for parametric densities. Let /i(x; 6) be 
a density on X that is parameterized by a parameter 0 in a 


parameter space 0. We extend /i(x; 6) to a (/-invariant density 
/ by using Thm. 2.1 obtaining: 


/( X; 0) = X^( X;6> )’ 


M 


M 


( 2 ) 


where /i,(x; 6) = h(GiX; 6). This density is of the form of a 
finite mixture of densities /ij(x; 6) of known parametric form 
where the mixture coefficients are all identical and equal to 
1 /M. Maximum likelihood (ML) estimation of the parameter 
0 from an i.i.d. sample {x,}” =1 from any (/-invariant density 
/ can now be performed using finite mixture model methods 
0 such as the Expectation-Maximization (EM) algorithm [2J 
or the restricted Boltzman machine (RBM) 0. 


III. ML within a Spherical Symmetry Group 

As in [Tj we specialize Thm. 2.1 to estimation of parameters 
for the case that the probability density is on a sphere and 
is invariant to actions in a spherical symmetry group. In 
Section [V] this will be applied to a crystallography example 
under spherical distribution likelihood models for the mean 
crystal orientation. In general, the measured and mean ori¬ 
entations can be represented by Euler angles [8], Rodrigues 
Vectors |9|, or Quaternions (TO). As in [ 11, we use the 
quaternion representation to enable orientations to be modeled 
by spherical distributions since the quaternion representation 
is a 4D vector on the 3-sphere S 3 , i.e. q = (qi, 92 , <Z 3 , 94 ) such 
that ||q|| = 1. 

Any of the aforementioned orientation representations have 
inherent ambiguity due to crystal symmetries. For example, if 
the crystal has cubic symmetry, its orientation is only uniquely 
defined up to a 24-fold set of proper rotations of the cube 
about its symmetry axes. These actions form a point sym¬ 
metry group, called 432, a sub-group of m3m. In quaternion 
space, since each orientation corresponds to two quaternions 
with different sign {q, — q}, these rotations reflections, and 
inversions can be represented as a spherical symmetry group Q 
of quaternionic matrices {Pi,..., Pat}, with sign symmetry 
such that P, = — Pi_M /2 VM/2 < i < M, where M = 48 
for cubic symmetry. 

Based on the symmetry group (/, we can define the distance 
between two quaternions under Q as: 


(qi t CI 2 ) = mjnarccos (q{Pq 2 ) (3) 

Two quaternions qi, q 2 are called symmetry-equivalent to 
each other if they are mapped to an equivalent orientation 
under (/, i.e. dg (qi • q 2 ) = 0. A fundamental zone (FZ), also 
called the fundamental domain, is a conic solid subset of the 
sphere that can be specified to disambiguate any particular 
orientation x. However, as will be seen in Sec.[V] reduction of 
the entire data sample {x, }'' =1 to a FZ destroys information 
necessary for maximum likelihood estimation: the entire Q- 
invariant density 0 must be used. In the following two 
subsections, we introduce two (/-invariant spherical distribu¬ 
tions: von Mises-Fisher and Watson distributions dj to model 
orientations in quaternion space. 







A. Hyperbolic Q-invariant von Mises-Fisher Distribution 

First we review the (/-invariant von Mises-Fisher distribu¬ 
tion presented in (I |. The von Mises-Fisher (VMF) distribution 
arises in directional statistics 03 as an analogue of the 
multivariate Gaussian distribution on the (p — 1)-dimensional 
sphere S'O’ -1 ) C R p , where p > 2. The VMF distribution 
is parameterized by the mean direction p, £ S (jl ~ ] ' and the 
concentration parameter n > 0 : 


4>(x; p, k ) = c p (k) exp (np T x), 


(4) 


where c v (n) = -rx and /„(•) is the modified 

Bessel function of the first kind of order p. Given an i.i.d. 
sample {x, }" = t from the VMF distribution, the ML estimator 
has the closed-form expressions CD 



l p! 2 


(n) 


where 7 = Ei=i x* and A p (u) = Ip/2 _^ v 

Let Q be a group of symmetric actions {Pi,...,P m} 
acting on the quaternionic representation of orientation on the 
3-dimensional sphere S 3 . We extend the VMF distribution (|4]i 


using the mixture representation in Thm 2.1 


/„(x;/x,k) 


1 M 

J] E 

( 6 ) 

1 M 

jj E Pm ^> k ) 

(7) 


m—1 


where in going from ([ 6 ]) to ((TJ) we used the inner product form 
fi ! x in (jdji and the symmetry of P m . The expression ([ 7 ]) for 
the extended VMF distribution is in the form of a finite mixture 
of standard VMF distributions on the same random variable x 
having different mean parameters p m = P m //, but having the 
same concentration parameter n. 

The finite mixture ([7|i for the (/-invariant density /„(x; p, k) 
is in a form for which an EM algorithm [|2j can be imple¬ 
mented to compute the ML estimates of p and n. Denot¬ 
ing the parameter pair as u = { p , k} , the EM algorithm 
generates a sequence {uE 7 } of estimates that monotonically 
increase the likelihood. These estimates are given by uj (k+l > = 
argmax^ [logL(oi; {x,, Sy})], where Si is a latent 

variable assigning x, to a particular mixture component in ([7]» 
and L(uj; {xi, Si}) is the likelihood function of uj given the 
complete data {x/,Si}” =1 . Specifically, 


where 77 TO = 
form: 

E-step: 

^i,m = 


£s|x,E lo g L ( w ;{ x i,Si})] (8) 

n M 

EE Cp(hv) + Av(P m /i.) Xi), 

i—l m—1 

P{ Si = m\xi,<jj). The EM algorithm takes the 


<ft(xj;P m /x,K) 


e {1,2,... ,M} . (9) 


M-step: 



(10) 

n M 


7 = E E G,m P ™ x i- 

(11) 


i—l m—1 


The hyberbolic Q -invariant von Mises-Fisher distribution is 
obtained by exploiting the sign symmetry in Q. In particular, 
(111 in the M-step can be re-written as: 

n M 

7 = E E 


i=l m—1 
n / M/2 

^ ] I y 1 r i,m P m x i — 

i= 1 \ m=l 
n M/2 


( 12 ) 


sinh (k(P mP) x^ 
tZ 1 EfE cosh («( p iF) T *i) 


= EE 


P T x 


where sinh, cosh are the hyperbolic sinusoidal functions. 
Equation (j9]l in E-step is simplified as: 


sinh (k(P m p) T Xj) 
cosl1 («(P ip) Tx i) 


m£ < 1,2,.... 


M 


(13) 

In Section [VI] we demonstrate the computational improve¬ 
ment of the hyperbolic form of the EM algorithm obtained by 
substituting l[T2}, <(T3]» into ([TOj, Q respectively. 


B. Q-invariant Watson Distribution 

As described at the beginning of this section, each orienta¬ 
tion corresponds to two quaternions with different sign, which 
is equivalent to an axis of the sphere. For axial data it is more 
natural to use the Watson distribution [12], which models the 
probability distribution of axially symmetric vectors on the 
(p—1)-dimensional unit sphere, i.e. ±x £ S p ~ 1 are equivalent. 
Similar to VMF, the distribution is parametrized by a mean 
direction ±p £ S p ~ 1 , and a concentration parameter k £ R 
that is no longer necessarily non-negative. Its probability 
density function is 

W p (x;p,K ) = T 1 p —-exp (n(p T x) 2 ), (14) 


where M is the Kummer confluent hypergeometric function 
defined in 03- According to ( [14} , the positive-negative pair 
of group actions {P m , — P m } contribute the same value in 
the density function. The set of the group action pairs Q' = 
{{P m ,-P m }}^i is the quotient group Q/ X, where 1 = 
{I p , —Ip} C G and Ip is the identity matrix of dimension 
p. Therefore, Q' is also a group and we can use Thm 2.1 to 


extend the Watson distribution to the mixture representation 
under Q'\ 


f w (x-p,K) 


1 

~M' 


M' 

^ ^ Wp(x; k), 

m—1 


(15) 











where M' = M/2. The ML estimates of //, and k can also be 
calculated by the EM algorithm. The E-step for the Watson 
mixture distribution is 


T i,m 


exp{K((P m /x ) r x,) 2 } 

E^' 1 exp{ K ((P iA x)^x i ) 2 } 


me {1,2,..., M'}. 


(16) 

For the M-step, we take a similar approach as GD as 
follows: 


^s|x,^[log L(u; {xj, Si})] 

=n^Kfj, T f/i-log ^M'M > (17 ' > 


where f = J E"=i Em=i is the scatter 

matrix of x 1; ..., x„. Let t\, ..., t p be the eigenvalues of T with 


f i > ... > tp , 


(18) 


and let ±ti,...,±t p be the corresponding unit eigenvectors. 
Since we want to find //, which maximizes GD such that 
[i 1 fj, = 1, the estimator of // for fixed k has the following 
form: 


(i = t 1; k > 0, 
jj, = t p , k < 0. 


(19) 


Similarly by fixing //, and setting to zero the derivative of 
with respect to re, we have: 


Consider the situation where the samples {x,}" = i follow 
a mixture of (/-invariant density functions. For the VMF 
distribution, the mixture density has the following form: 


M 


3 „(x; {n c , k c , a c }) = a c ( K c) 


C— 1 


\m—l 


( 21 ) 

where C is the number of clusters assumed to be fixed a priori, 
/x c) re c are the parameters for the c-th cluster and a c are the 
mixing coefficients where E^L i a c — 1 an< l a c > 0 for all c. 
The parameters of (21 1 can be estimated by the EM algorithm: 
E-step: 


PmMc) «c) 


Eh=i a h Ei=i 0(xi; P iHh, Kh.) 


( 22 ) 


M-step: 


OL C 


7 C 


M 


— ||^ J , «c — A 


-1 

p 


2—1 m—1 
n M 


^ 1 ^ ] r i,c,mP m x !) 


2—1 m=l 



, (23) 
(24) 


where r, jC . m is the probability of sample x, belonging to the 
c-th cluster and the m-th symmetric component. 

For the Watson distribution, the mixture of (/-invariant 
Watson density is 


Y p {k) = 


ni^*)_Er=i£:=t^(M T PM 


M' 




/Y' n r (tt T P T x- 

-1 ( Z_/i=l Z^m=l 1 h rn \t J ‘ r m-^i 


n 

TnT „ \2 ' 


=>re = Y1 


( 20 ) 

The final estimates of // and re are obtained by checking both 
cases (re > 0, re < 0) and choosing the one which is consistent 
for <fH 


IV. Clustering with a Spherical Symmetry Group 

In this section we extend the parameter estimation problem 
to the situation where there are multiple group-invariant distri¬ 
butions with different parameters that govern the samples. This 
problem arises, for example, in poly-crystaline materials when 
estimating the mean crystal orientation over a region contain¬ 
ing more than one grain (perhaps undetected). This problem 
can be solved by first applying some standard clustering meth¬ 
ods, e.g. K-means 03’ and then estimating the parameters 
for each cluster. However, clustering methods based on the 
distance relation between the samples are complicated by the 
presence of spherical symmetry because it is necessary to dis¬ 
tinguish modes that are due only to symmetry from those that 
distinguish different clusters. Therefore, we propose a model- 
based clustering algorithm which accommodates symmetry to 
handle this problem. 


C 




, (x; {/r c , k c , a c }) — ^ a c j ^ ' ^^ / Wp(x.',P m /i c , 


c= 1 


t m =1 


The E-step is similar to (22 1 with d replaced by W p 
function. The M-step can be computed with a similar approach 
as in Section |III-B| with the following modifications: 


n M r 


Tn = 


-zL r, i,c,m(P^X i xf p m ) J 

na c Lz ' 


k, = Y„ 


/u" r (u T P T x-l 2 ' 

-1 / Z^i=lZjm=l'wlrc r m A V 


(26) 

(27) 


1 UAH 

where a c = E»=i E m =i G,c,m- 


A. Multi-modality Tests on Q-invariant Spherical Distributions 

Given sample set (x, }" = i on S' 1 * -1 , the objective is to 
determine whether the n samples are drawn from one single 
distribution or a mixture of C distributions. For polycrystalline 
materials, the result of this determination can be used to 
discover undetected grains within a region. We propose to 
use a multi-modal hypothesis test based on the (/-invariant 
distributions to solve this problem. The two hypotheses are 
//(,: The samples are from a single (/-invariant distribution 
/(x; {/x, re}); and H \: The samples are from a mixture of 














C distributions g(x; {p c , k c . a,,}/+, )• The Generalized Like¬ 
lihood Ratio Test (GLRT) |15 | has the following form: 

{^t e ,K e ,a c }g =1 G0i ff({ x i}£= li {Mci Q!c}^=l) 


A, 


max 


GLR = 


max {MjK}e e 0 /({xi}-U; {/+ «}) 




(28) 

where 0 O , 0 X are the parameter spaces for the two hypotheses. 
The / and g functions for VMF and Watson distributions 
are defined in ([7]), ( p~5j ) and (21 1 , ( |25} respectively and the 
test statistic A glr can be calculated by the proposed EM 
algorithm. According to Wilks’s theorem [16] as n approaches 
oo, the test statistic 2 log A glr will be asymptotically x 2 - 
distributed with degrees of freedom equal to (jp + 1 )(C — 1), 
which is the difference in dimensionality of 0o and 0i. 
Therefore, the threshold ?/ in ( |28| can be determined by a 
given significance level a. 

V. Application to Crystallographic Orientation 

Crystal orientation and the grain distribution in polycrys¬ 
talline materials determine the mechanical properties of the 
material, such as, stiffness, elasticity, and deformability. Lo¬ 
cating the grain regions and estimating their orientation and 
dispersion play an essential role in detecting anomalies and 
vulnerable parts of materials. 

Electron backscatter diffraction (EBSD) microscopy ac¬ 
quires crystal orientation at multiple locations within a grain 
by capturing the Kikuchi diffraction patterns of the backscatter 
electrons 03- A Kikuchi pattern can be translated to crystal 
orientation through Hough Transformation analysis 1181 or 


Dictionary-Based indexing GD- The process of assigning 
mean orientation values to each grain is known as indexing. 
Crystal forms possess point symmetries, e.g. triclinic, tetrag¬ 
onal, or cubic, leading to a probability density of measured 
orientations that is invariant over an associated spherical 
symmetry group Q. Therefore, when the type of material has 
known symmetries, e.g., cubic-type symmetry for nickel or 
gold, the ^-invariant VMF and Watson models introduced in 


Section III can be applied to estimate the mean orientation 


p CJ and the concentration n g associated with each grain. 
Furthermore, the clustering method along with the multi¬ 
sample hypothesis test in Section [TV] can be used to detect 
the underlying grains within a region. 


A. Simulation of Crystallographic Orientation 

To simulate the crystallographic orientations, we first draw 
random samples from VMF and Watson distributions with p = 
4. The random variable x in a spherical distribution can be 
decomposed 0 = 

x = f/it + \/l - tfSf (x), (29) 

where t = p T x and (x) = (I p — pp T )x/\\(I p — pp T )x\\. 
Let /(x; /it) be the p.d.f. of the distribution where p is 
the mean direction. According to the normal-tangent decom¬ 
position property, for any rotationally symmetric distribu¬ 
tion, Sfj.(x.) is uniformly distributed on , the (p — 2)- 


dimensional sphere normal to p, and the density of t = x T p 
is given by: 

t /c/(f)(l-f 2 ) (p ~ 3)/2 . (30) 

For VMF distribution, substituting ( [29} into 0 and combin¬ 


ing with (30 1 , we have the density of the tangent component 
t as: 


fv(t) = C„exp{/ct}(l - t 2 ) (p 3)/2 

C v v 2 / WlW r 


r '^ 


-l 


(31) 


Similarly, the density of the tangent component of Watson 
distribution is: 


fwit) = C w exp{Kt 2 }(l -f 2 ) (p 3)/2 


C w = 


m 


(32) 


r( ^)r(I)M(i,I/O- 


Random samples from the density functions ( fTi~| and (32 1 can 
be easily generated by rejection sampling. 

The generated quaternions from VMF and Watson 
distributions are then mapped into the Fundamental Zone 
(FZ) with the symmetric group actions to simulate the wrap¬ 
around problem we observe in real data, i.e. observations 
are restricted to a single FZ. For cubic symmetry, the FZ in 
quaternion space is defined in the following set of equations: 


1 92/91 1 < y / 2 - 1 
1 93 / 51 1 < \/2- 1 

|94/gi| < \/2 - 1 

I92/91 + 93M +94/91 1 < 1 
I92/91 - 93/91 +94/91 1 < 1 
1 92/91 + 93/91 - 94/91 1 < 1 
1 92/91 - 93/91 - 94/91 1 < 1 

where q, is the ?-th component of quaternion q. 


I92/91 - 93/91 1 < x /2 
I92/91 + 93/91! < x /2 
I92/91 - 94/91! < x /2 
I92/91 + 94/91 1 < x /2 
I93/91 - 94/91 1 < \/2 
I93/91 + 94/91 1 < x /2 


(33) 


VI. Experimental Results 

A. Q-invariant EM-ML Parameter Estimation on Simulated 
Data 

Sets of n i.i.d. samples were simulated from the VMF or 


Watson distributions using the method described in Sec V-A 
with given p = p a ,n = k 0 for the to3to point symme¬ 
try group associated with the symmetries of cubic crystal 
lattice planes. The number of samples for each simulation 
was set to n = 1000 and n 0 was swept from 1 to 100 
while, for each simulation run, fi a was selected uniformly 
at random. The experiment was repeated 100 times and the 
average values of k and the inner product p p a are shown 
in Fig. [I] and [2] In the figures we compare performance 
for the following methods: (1) the naive ML estimator for 
the standard VMF or Watson model that does not account 
for the point group structure (labeled ”ML Estimator”). (2) 
Mapping each of the samples x, toward a reference direction 
x r (randomly selected from {x;}" =1 ), i.e. x, 1 —»• P m x t , 
where P m = argminp e p arccos (x/Px), to remove possible 


















(a) VMF Simulated Data 



(b) Watson Simulated Data 




(a) VMF Simulated Data 


(b) Watson Simulated Data 


Fig. 1. Mean orientation estimator comparisons for (/-invariant densities. 
Shown is the average inner product /i ^ ft of four estimators fi when /x G is 
the true mean orientation as a function of the true concentration parameter 
k 0 for the data simulated from VMF (Fig. |l(a)) and from Watson (Fig. Eg} 
distribution. The naive estimator (”ML Estimator” in blue line) does not attain 
perfect estimation (inner product = 1) for any k 0 since it does not account 
for the spherical symmetry group structure. The modified ML (green dashed 
line) achieves perfect estimation as k q becomes large. The proposed EM-ML 
methods (’’EM-VMF”, ’’EM-Watson”) achieve perfect estimation much faster 
than the other methods even under model mismatch (EM-VMF for Watson 
simulated data and vice versa). 


ambiguity. Then performing ML for the standard VMF or 
Watson distribution (labeled ’’Modified ML”). (3) Applying 
our proposed EM algorithm directly to the n samples using 
the mixture of VMF distribution (|9|)-(fTT]) (labeled ’’EM-VMF”) 
(4) Applying our proposed EM algorithm to the mixture of 
Watson distribution (p~6]>-([20]> (labeled ’’EM-Watson”). 


Figure [I] shows the inner product values fi[ fi. The proposed 
EM-VMF and EM-Watson estimators have similar perfor¬ 
mance in that they achieve perfect recovery of the mean 
orientation ( fi,]~ fi = 1) much faster than the other methods as 
the concentration parameter k 0 increases (lower dispersion of 
the samples about the mean) no matter whether the data is gen¬ 
erated from VMF (Fig. 1 1 (a) i or Watson distribution (Fig. 1 (b) i, 
indicating the robustness of the proposed approaches under 
model mismatch. Notice that when n 0 is small (k 0 < 20 
for VMF data and k 0 < 10 for Watson data), none of the 
methods can accurately estimate the mean orientation. The 
reason is that when k 0 is small the samples become nearly 
uniformly distributed over the sphere. The threshold k 0 value 
at which performance starts to degrade depends on the choice 
of point symmetry group and the distribution used to simulate 
the data. In Fig. [2] it is seen that the biases of the proposed 
EM-VMF 11 ] and EM-Watson n estimators are significantly 
lower than that of the other methods compared. While the 
modified ML performs better than the naive ML estimator, its 
bias is significantly worse than the proposed EM-VMF and 
EM-Watson approaches. 


Figure [3] shows the computation time of the estimation 
algorithms presented in Fig. [I] and Fig. [2] The computation 
time for all methods decreases as k 0 becomes larger. When 
k 0 is small (k 0 < 20 for VMF data and k 0 < 10 for Watson 
data), because the samples are almost uniformly distributed 
around the sphere, it is difficult for the EM algorithms to 
converge to the optimal solution and they therefore require 
maximum number of iterations to stop, forming the plateaus 


Fig. 2. Concentration parameter estimator bias as a function of the true 
concentration k, 0 for data simulated from VMF (Fig. |2(a)) Q and from 
Watson (Fig. [2(b)} distributions. The bias of the naive ML (blue solid line) is 
large over the full range of k 0 . The modified ML (green dashed line) estimates 
k, more accurately when k 0 is small. Our proposed EM-VMF and EM-Watson 
estimators (black dotted line and magenta dashed line) have lower bias than 
the other estimators. 
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Fig. 3. Computation time for calculating the result in Fig. |T] and Fig. [2] EM- 
Watson (magenta dashed line) has less computation time than EM-VMF (black 
dotted line) because it uses only half of the symmetry operators. EM-VMF- 
Hyper (cyan circle line) which uses the hyperbolic sinusoidal simplification 
of EM-VMF reduces the computation time by more than a factor of 2. 


in Fig. [3] Notice that EM-Watson requires less time than EM- 
VMF even though it has more complicated E and M-steps. 
The reason is that EM-Watson uses only half of the symmetry 
operators, which corresponds to the size of the quotient group 
Q /X as described in Section III-B By applying the hyperbolic 
sinusoidal simplification in Section [III-A (labeled ’’EM-VMF- 
Hyper”), we can further reduce the computation time by more 
than a factor of 2 compared to the original EM-VMF. 


B. Q-invariant Clustering on Simulated Data 

In this section, we demonstrate the performance of our pro¬ 
posed EM approaches for clustering. Sets of n i.i.d. samples 
were simulated from the VMF or Watson distributions with 
k = k 0 and one of two mean directions (//,, fi 2 ) to generate 
two clusters of samples. The spherical symmetry group is 
m3m as before. The number of samples for each set was 
set to n = 1000 and k 0 was swept from 1 to 100 while, 
for each set, . fi 2 was selected uniformly at random. The 
experiment was repeated 100 times and the average values of 
the inner product (/if Hi + fh M 2 V 2 are shown in Fig. [T] In 
the figure we compare performances of the following methods: 
(1) Cluster the samples by standard K-means algorithm with 
the distance defined by the arc-cosine of the inner product 






















































and then use the naive ML within each cluster to estimate the 
mean directions (labeled ”K-means”). (2) Cluster the samples 
by K-means with the distance defined as 0 and then use 
the aforementioned modified ML estimator (labeled ’’Modified 
K-means”). (3) Apply our proposed multi-cluster EM-VMF 
algorithm to the n samples directly (|22)-(24i (labeled ’’EM- 


VMF”) (4) Apply our multi-cluster EM-Watson algorithm to 
the n samples directly @-(127) (labeled ’’EM-Watson”). 

Figure shows the average inner product values + 

As"A *2 )/2 from the mean direction estimation. The proposed 
EM-VMF and EM-Watson are able to correctly cluster the 
samples and achieve perfect recovery of the two mean orien¬ 
tations much faster than the other K-means approaches. Notice 
that the region where all the methods fail is larger than the 
single cluster case since multiple clusters increase the difficulty 
of parameter estimation. Again, no matter whether the samples 
are simulated from VMF or Watson distribution, our proposed 
approaches perform equally well under both cases. 

To further test the ability to detect multiple clusters given 
a set of samples, we generate 1000 sets of samples. Each set 
has 1000 samples and is assigned randomly to label 0 or 1. If 
the set is labeled 0, the samples are generated from a single 
distribution; If the set is labeled 1, then the samples in the set 
are randomly generated from two distributions with different 
means. The GLRT is used with the four aforementioned clus¬ 
tering methods to test whether the samples in each set are uni- 
modal or multi-modal. The Receiver Operating Characteristic 
(ROC) curves of the test results are shown in Fig. [5] The 
naive K-means with ML estimator which does not consider 
the symmetry group actions fails to distinguish whether the 
multiple modes are from actual multiple distributions or due 
to the wrap-around effect from the fundamental zone mapping. 
Therefore, this approach tends to over-estimate the goodness 
of fit of the Hi model for true negative cases and under¬ 
estimate it for true positive cases, resulting in a result that 
is even worse than random guessing. The modified K-means 
performs better than K-means but worse than our proposed 
EM-VMF and EM-Watson algorithms. 


C. EM-ML orientation estimation for IN100 Nickel Sample 

We next illustrate the proposed EM-VMF and EM-Watson 
orientation estimators on a real IN 100 sample acquired from 
US Air Force Research Laboratory (AFRL) |T9] |. The IN 100 
sample is a polycrystalline Ni superalloy which has cubic 
symmetry in the rn'-irn point symmetry group. EBSD orienta¬ 
tion measurements were acquired on a 512 x 384 pixel grid, 
corresponding to spatial resolution of 297.7 nm. The Kikuchi 
diffraction patterns were recorded on a 80 x 60 photosensitive 
detector for each of the pixels. 

Figure [6] (a) shows a 200 x 200 sub-region of the full EBSD 
sample where the orientations are shown in the inverse pole 
figure (IPF) coloring obtained from the OEM EBSD imaging 
software and (b) is the back-scattered electron (BSE) image. 
Note that the OEM-estimated orientations in some grain 
regions of the IPF image are very inhomogeneous, having 
a mottled appearance, which is likely due to a fundamental 
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Fig. 4. Mean orientation estimator comparisons for samples generated from 
two different means. Shown is the average inner product (p,J Mi + Ai M g)/2 
of four methods when M 11 M 2 are the true mean orientations as a function 
of the true concentration parameter k 0 for the data simulated from VMF 
(Fig. |4(aj) and from Watson (Fig. |4(bj) distributions. The K-means with naive 
estimator (”K-means” in blue line) does not attain perfect estimation for any 
k 0 . A modified K-means with ML estimator (’’modified K-means” in green 
dashed line) achieve perfect estimation as k 0 becomes large. The proposed 
EM-VMF and EM-Watson methods ("EM-VMF”, ”EM-Watson”) achieves 
perfect estimation much faster than the other methods. 
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Fig. 5. ROC curve for detecting bi-modally distributed samples. The samples 
are uni-modal or bi-modal distributed from VMF (Fig. |5(a)| or Watson 
( |5(b)| distributions with k q = 50. The naive K-means with ML estimator 
cannot cluster the samples well and estimate the mean directions accurately, 
resulting in poor detection which is even worse than random guessing. The 
modified K-means (green dashed line) performs better than K-means but 
is still unsatisfactory. Our proposed EM-VMF (black dots) and EM-Watson 
(magenta dashed line) methods have very good performance in this detection 
task. 


zone wrap-around problem. As an alternative, we apply a 
combination of the proposed EM estimators (EM-VMF or EM- 
Watson) and the GLRT ( [28) with C = 2 and significance 
level a = 0.05 to detect multi-modal distributions within each 
OEM-segmented region. Figure [6] (c)(e) show the estimates 
of the mean orientations of the regions/sub-regions, where the 
sub-regions surrounded by white boundaries indicate those that 
have been detected as deviating from the distribution of the 
majority of samples from the same region. The multi-modally 
distributed regions may be due to undetected grains, inaccurate 
segmentation, or noisy orientation observations. To distinguish 
the latter situations from the first in which the region really 
consists of two grains, the misalignment/noise test introduced 
in |[20) can be used. Figures [6] (d)(f) show the estimated con¬ 
centration parameter k for the regions/sub-regions. Note that 
the estimated k are large for most of the regions/sub-regions 
because those regions which have multi-modally distributed 


























samples are detected and their concentration parameters are 
estimated separately for each sub-region. 
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Fig. 6. A 200 x 200 sub-region of the IN 100 sample, (a) is the IPF 
image for the Euler angles extracted from EBSD by OEM imaging software. 
IPF coloring in some grains is not homogeneous, likely due to symmetry 
ambiguity, (b) is the BSE image of the sample. (c)(e) show the estimates 
of the mean orientations of the regions/sub-regions using a combination 
of the proposed EM estimators, EM-VMF and EM-Watson respectively, 
and the GLRT i |28) to detect multi-modal distributions within each OEM- 
segmented region. The sub-regions surrounded by white boundaries indicate 
those that have been detected as deviating from the distribution of the majority 
of samples from the same region. (d)(f) show the estimated concentration 
parameter k for the regions/sub-regions. Note that the estimated k, are large 
for most of the regions/sub-regions because those regions which have multi- 
modally distributed samples are detected and their concentration parameters 
are estimated separately for each sub-region. 


VII. Conclusion 

A hyperbolic Q -invariant von Mises-Fisher distribution was 
shown to be equivalent to the distribution proposed in |lj. 
The advantage of the hyperbolic form is parameter estimation 
can be performed with substantially fewer computations. A 
different group invariant orientation distribution was intro¬ 
duced, called the ^-invariant Watson distribution, and an EM 
algorithm was presented that iteratively estimates its orienta¬ 
tion and concentration parameters. We introduced multi-modal 
generalizations of these ^-invariant distributions using mixture 
models and showed that these can be used to effectively cluster 


populations of orientations that have spherical symmetry group 
invariances. The mixture of VMF and Watson models were 
applied to the problem of estimation of mean grain orientation 
parameters in polycrystalline materials whose orientations lie 
in the m3m point symmetry group. Application of the finite 
mixture representation to other types of groups would be 
worthwhile future work. 
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